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Introduction 

The  overall  technical  objective  of  the  Phase  I  work  was  to  demonstrate  the  feasibility  of  a  compu¬ 
tational  virtual  laboratory  for  simulating  high  voltage  effects,  herein  called  VLAB.  Models  in  sup¬ 
port  of  this  objective  were  identified  and  defined  in  Phase  I,  with  some  components  implemented 
in  the  one-dimensional  object  oriented  code,  OOPD1,  and  others  implemented  directly  in  the  2D 
code  XOOPIC.  OOPD1  serves  as  a  training  vehicle  for  new  students  and  a  platform  for  rapid  pro¬ 
totyping  of  new  algorithms.  The  simplicity  of  OOPD1  ensures  that  the  algorithm  or  model 
remains  the  focus,  rather  than  the  mathematical  and  programming  details  of  multi-dimensional 
implementation.  In  Phase  II,  all  final  implementation  will  be  in  the  two-dimensional  VLAB  code, 
based  on  the  well-tested  XOOPIC  particle-in-cell  code. 

The  results  of  Phase  I  of  this  work  are  summarized  in  this  section.  In  Task  1,  the  work  related  to 
modeling  the  formation  of  a  moving  plasma  cathode  interface  is  described.  The  Task  2  work 
related  to  scattering  of  energetic  electrons  is  discussed.  Work  related  to  high  voltage  breakdown 
of  insulator  surfaces  is  discussed  in  Task  3.  Models  for  intense  heat  fluxes  to  surfaces  is  described 
in  Task  4,  and  models  for  generation  of  x-rays  are  discussed  in  Task  5.  Work  related  to  detailed 
structural  visualization  of  high  voltage  microwave  devices  is  discussed  in  Task  6. 

Task  1:  Moving  Plasma  Cathode 

Plasma  formation  near  electrode  surfaces  can  result  when  desorbed  impurities  or  sputtered  neu¬ 
trals  are  ionized  by  energetic  electrons  or  x-rays.  Even  in  the  presence  of  strong  fields,  the  space 
charge  of  the  plasma  can  lead  to  formation  of  a  virtual  cathode.  Temporal  variations  in  the  fields 
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or  plasma  density  can  lead  to  translation  of  the  point  at  which  the  field  is  zero,  leading  to  the  anal¬ 
ogy  of  a  moving  Child-Langmuir  space  charge  limited  emitter  in  the  presence  of  fields  that  extract 
electrons  from  the  plasma. 

A  suite  of  models  was  developed  for  emission  from  a  moving  plasma  cathode  surface  during 
Phase  I  of  the  research.  The  following  models  are  required  for  study  of  emission  from  plasma 
cathodes.  An  effective  field  enhancement  model  was  developed  to  incorporate  geometric  surface 
effects  at  the  sub-grid  level,  as  well  as  to  include  geometric  field  enhancement  in  models  with 
lower  spatial  dimensionality.  The  transition  from  Fowler-Nordheim  field  emission  current  to 
space-charge  limited  current  was  also  described.  We  also  investigated  hybrid  models  for  repre¬ 
senting  more  dense  plasmas,  which  may  be  formed  near  emitting  surfaces,  particularly  particle- 
fluid  hybrids  capable  of  retaining  important  kinetic  effects.  The  standard  Monte  Carlo  collision 
model  was  extended  to  the  relativistic  regime.  Sputtering  and  desorption  models  for  particle 
impact  at  surfaces  were  also  described. 

Geometric  Field  Enhancement  of  Emission 

A  field  enhancement  factor  is  often  required  at  the  macroscopic  scale  in  order  to  explain  field 
emission  current  densities  due  to  geometric  features.  The  local  field  enhancement  factor  p  is  often 
introduced  in  the  Fowler-Nordheim  equation  to  represent  the  geometrical  effects  at  the  surface  of 
the  cathode,  where  P(s)  =  E(s)/E0 .  Local  variation  of  P  determines  the  local  surface  electric 

field,  resulting  in  local  dependence  of  injection  current  by  the  Fowler-Nordheim  law.  In  computa¬ 
tional  models,  it  is  impractical  to  determine  the  time-dependent  local  surface  field  for  each  time 
step  on  a  microscopic  space  scale.  Effective  P  is  introduced  in  order  to  facilitate  study  of  the  emis¬ 
sion  properties  at  a  macroscopic  scale.  Microscopic  (subgrid)  local  p  is  calculated  only  at  the  ini¬ 
tial  time  step,  and  then  the  effective  enhancement,  Pe^,  can  be  recomputed  under  different  surface 
electrical  fields  through  this  model.  The  model  allows  reduction  of  dimensionality  as  well  as  the 
ability  to  include  subgrid  effects.  The  model  is  demonstrated  on  fundamental  cases  and  compared 
to  a  calculation  with  a  mesh  fine  enough  to  resolve  the  geometric  features. 

The  basic  field  emission  effect  involves  reduction  of  the  tunneling  barrier  for  electrons  within  a 
surface,  described  by  the  well-known  Fowler-Nordheim  equation,  modified  here  to  include  the 
local  field  enhancement  factor  [1]: 
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where  E(s)  is  the  normal  component  of  electric  field  at  a  ID,  2D,  or  3D  surface  s ,  <j)w  is  the 

work  function  of  the  metal,  A  and  B  are  empirical  constants,  y  =  3.79x10  (P(s)2?0(j))  , 

and  p(s)  =  E(s)/Eg(s)  is  the  ratio  of  the  local  surface  electric  field  to  the  applied  field. 

\)(y)=0.95-y2  and  t(y)  «  1.1  are  functions  that  arise  due  to  the  inclusion  of  image  charge  effects.  It 
has  been  noted  that  this  field  enhancement  may  come  not  only  from  geometric  effects  but  also 
from  impurities  [2,3]. 

The  methodology  for  incorporating  pe^ris  described  by  the  following.  Starting  from  the  surface 
geometry,  the  local  electric  fields  due  to  an  applied  field  are  solved  at  t= 0  to  obtain  the  local  field 
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enhancement  (3  at  each  point  along  the  surface.  Next,  the  effective  current  for  a  given  applied  field 
is  found  and  then  Pe^is  determined  over  the  emitting  surface. 

Effective  field  enhancement  p^is  defined  for  most  the  general  case  by  integrating  the  current 
over  the  emitting  surface  numerically, 
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where  As  is  the  emitting  surface  area.  In  our  implementation,  we  subdivide  the  cathode  with  many 
small  intervals  and  use  Simpson’s  3/8  rule  to  integrate  (2)  in  each  interval,  then  use  the  bisection 
method  to  invert  (3)  for  p ejp  For  the  case  that  E0(s)  is  a  constant  over  the  surface  of  cathode,  we 
can  simply  set 
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Without  the  consideration  of  AE,  the  effective  P  will  include  all  the  information  for  the  surface  of 
the  cathode.  This  can  be  applied  to  reduce  the  dimensionality  required  to  model  higher  dimension 
geometry,  as  well  as  allow  for  gridding  which  does  not  resolve  the  underlying  geometry.  The  pro¬ 
cess  is  illustrated  by  example  in  Fig.  1  for  an  analytic  local  field  enhancement  factor. 


Fowler-Nordheim  Transition  to  Space  Charge  Limited  Emission 

The  Fowler-Nordheim  law  gives  the  current  density  extracted  from  a  surface  under  strong  fields 
by  treating  the  emission  of  electrons  from  a  metal-vacuum  interface  in  the  presence  of  an  electric 
field  normal  to  the  surface  as  a  quantum  mechanical  tunneling  process.  Child’s  law  predicts  the 
maximum  transmitted  current  density  by  considering  the  space  charge  effect.  At  low  fields,  the 
transmitted  current  is  limited  by  Fowler-Nordheim’s  law,  and  at  high  fields  it  is  limited  by  Child’s 
law.  This  work  analyzes  the  transition  of  the  transmitted  current  density  from  the  Fowler-Nord¬ 
heim  law  to  the  Child’s  law  space  charge  limit  using  the  one-dimensional  OOPD1  code.  Also 
studied  is  the  response  of  the  OOPD1  emission  model  to  strong  electric  fields  near  the  transition 
point.  Methods  for  increasing  the  current  density  to  approach  the  space  charge  limit  by  modifying 
the  device  design  are  also  described. 

In  Fig.  2,  the  Fowler-Nordheim  current  is  compared  to  the  space  charge  limited  current  for  P=1 

and  p=10  in  a  fixed  10'6  m  gap  for  varying  applied  field.  Here  applied  field  is  defined  as  the 
applied  gap  voltage  divided  by  the  gap  width.  In  this  model,  the  p  is  held  fixed,  when  in  fact  it 
should  appear  as  p^and  vary  with  the  applied  field  according  to  the  preceding  section.  Note  that 
in  these  coordinates,  the  space  charge  limiting  current  is  independent  of  p,  and  depends  only  on 
the  applied  field  (with  gap  held  fixed,  and  holding  emission  velocity  fixed  and  small  compared  to 
the  applied  voltage).  From  these  results,  it  requires  P  >  10  to  approach  the  space  charge  limited 
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Figure  1.  Example  of  effective  beta: 
(a)  local  dependence  of  P,  (b) 
averaged  Fowler-Nordheim  current 
density,  and  (c)  fitted  line  between 
Pe^-and  surface  electric  field 
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current,  as  the  space  charge  otherwise  suppresses  the  surface  field,  and  neither  the  Fowler-Nord- 
heim  current  nor  the  space  charge  limited  current  can  be  reached  even  at  high  applied  field 
strength. 

Hybrid  Models 

One  significant  computational  limitation  of  modeling  the  formation  of  a  plasma  cathode  is  that 
the  densities  can  become  significant,  making  a  particle  approach  computationally  intensive.  How¬ 
ever,  in  order  to  capture  non-equilibrium  and  kinetic  behaviors,  a  particle  method  is  preferred.  In 
order  to  make  PIC  codes  more  efficient  for  high  density  distributions,  which  contain  a  significant 
thermal  part  of  the  distribution,  a  hybrid  particle-fluid  model  is  considered.  In  Phase  I,  a  number 
of  candidate  methods  were  identified,  as  discussed  below. 

Fluid  behavior  is  determined  by  the  Navier-Stokes  equations: 

+  V  •  (nv)  =  0  (EQ5) 

mn(^t  +  (\  ■  V)\j  =  qne(E  +  v  x  B)-mnxvm-  Vp  (EQ6) 


^)+V^(pv)+^v.v  +  v-Q 


(EQ7) 


Relativistic  forms  may  be  used  if 
necessary. 

The  ideal  gas  equation  is  used  as  an 
equation  of  state: 

p  =  nT  (EQ  8) 

The  temperature,  T,  is  in  eV. 

First,  consider  electrostatic  hybrid 
models.  One  simple  model,  often 
called  a  Boltzmann-PIC  hybrid  [4], 
represents  part  of  the  electron  distri¬ 
bution  with  inertialess  fluid  electrons 
obeying  Boltzmann’s  relation. 
Neglecting  the  left  side  of  Eq.  (6), 
and  assuming  driftless  species  with 
pressure  balancing  the  electric  field 
for  the  isothermal  case, 

0  =  nq'V®  +  V(neTr)  (EQ9) 

n(x)  =  «0exp(-^O(x)/T)  (EQ  10) 


Electric  Field  (V/m) 

Figure  2.  Transition  from  Fowler-Nordheim  to 
space  charge  limited  current  for  a  fixed  gap  of  10'6 
m,  for  p  =1  and  P  =10. 
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M0  is  the  density  where  0  =  0.  Poisson’s  equation  becomes  a  nonlinear  expression  for  O.  For  a  sin¬ 
gle  Boltzmann  species: 

V20  =  [pM/(x)  +  ^e»0exp(-gO)(x)/T)]  (EQ  11) 

£o 

perf(x)  is  the  charge  density  at  x  from  all  sources  other  than  the  Boltzmann  species.  Eq.  (11)  may 
be  solved  with  a  nonlinear  numerical  solver. 

The  next  level  of  complexity  of  the  fluid  model  is  the  drift-diffusion  approximation.  The  drift-dif¬ 
fusion  approximation  also  includes  the  same  two  right-hand  side  terms  from  the  momentum  equa¬ 
tion  (Eq.  (6))  as  the  Boltzmann  approximation  and  also  includes  the  collisional  term,  mn\vm : 
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q  E 


eT  Vn 
m\  n 
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(EQ  12) 


This  assumes  vm  is  known  or  may  be  calculated.  Assuming  T  is  known,  there  are  now  two  vari¬ 
ables  («,v),  as  opposed  to  one  for  the  Boltzmann  approximation.  The  system  may  be  closed  by 
using  the  continuity  equation  (Eq.  (5)).  The  charge  density  is  added  to  the  Poisson  equation,  as  in 
the  Boltzmann  approximation. 

Sommerer  and  Kushner  [5]  simulate  RF  discharges,  including  chemical  reactions  in  gas  mixtures 
(He,  N2,  02,  CF4,  SiH4,  NH3)  using  a  ID  hybrid  model  consisting  of  a  Poisson  solver  and  drift- 
diffusion  fluid  model  (including  chemistry  terms),  an  electron  Monte  Carlo  simulation,  and  a  neu¬ 
tral  chemistry  and  transport  model.  Ion  mobilities  are  from  tabulated  values  and  electron  mobili¬ 
ties  are  generated  as  a  curve  fit  to  the  collision  frequencies  calculated  from  the  previous  Monte 
Carlo  iteration.  The  diffusion  coefficient  is  determined  using  the  Einstein  relation,  D  =  \ikT.  Aver¬ 
aging  is  done  over  several  RF  cycles,  keeping  the  phase  of  the  dependent  variables,  and  the  elec¬ 
tron  Monte  Carlo  model,  the  fluid  model,  and  the  chemistry  model  are  iterated  over  until  the 
system  converges  to  an  oscillatory  system.  Results  are  compared  with  experimental  data. 
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Since  Coulomb  collisions  are  impor¬ 
tant  for  high  density  plasmas,  it  is 
problematic  that  the  standard  PIC 
method  under-represents  the  Coulomb 
force  for  particles  within  a  cell,  as 
shown  in  Fig.  3.  A  new  technique 
developed  partially  under  Phase  I  of 
this  program  is  the  treecode  method 
[6],  in  which  particles  interact  directly 
via  Coulomb’s  law  at  short  range,  but 
interact  as  clusters  of  particles  at  long 
range. 

The  homogenous  solution  is  obtained 
via  a  boundary  integral  method.  This 
allows  the  method  to  include  both  the 


short  and  long  range  Coulomb  colli¬ 
sion  operators,  while  still  resulting  in 
an  N  log  N  method  for  N  particles, 
comparable  to  the  PIC  scheme. 


Figure  3.  Force  between  two  PIC  particles  separated 
by  a  distance  r  using  linear  weighting,  compared  to 
Coulomb’s  law  in  ID,  2D,  and  3D 


Because  the  method  is  gridless,  it  offers  an  efficient  scheme  for  modeling  non-conformal  bound¬ 


aries.  The  method  requires  spatial  sorting  of  particles.  Many  details  have  not  yet  been  resolved  for 


the  treecode  method,  such  as  the  consequence  for  noise  properties  resulting  from  the  short  range 


Coulomb  interaction. 


A  comparison  of  the  treecode  scheme  with  the  classical  PIC  model  is  shown  in  Fig.  4  for  a  one¬ 
dimensional  virtual  cathode  oscillation.  Note  that  the  PIC  method  results  in  a  sharper  oscillation, 
because  the  short  range  Coulomb  collisions  in  the  treecode  result  in  lower  peak  densities  and  con¬ 
sequently  a  weaker  interaction  at  the  peak.  The  higher  peak  densities  at  the  turning  point  in  the 
PIC  code  result  in  exaggeration  of  the  current  modulation,  which  can  be  seen  as  gaps  in  the  veloc¬ 
ity-position  phase  space.  A  direct  application  of  Coulomb’s  law  agrees  well  with  the  treecode 
result,  while  the  PIC  result  converges  to  the  treecode  result  for  A  x— >0. 

Next,  we  consider  electromagnetic  models.  The  drift-diffusion  approximation  may  also  be  formu¬ 
lated  with  the  assumption  of  a  magnetic  field  for  electromagnetic  plasmas.  In  this  case,  there  are 
two  diffusion  and  mobility  terms  -  those  along  field  lines  and  those  across  field  lines.  The  quanti¬ 
ties  along  the  magnetic  field  are  the  same  as  the  electrostatic  case.  As  with  the  electrostatic  case, 
the  momentum  collision  frequency  must  be  known  or  calculable. 

The  mobility  and  diffusion  coefficients  for  crossing  magnetic  field  lines  are  related  to  electrostatic 
quantities  [7]: 


(EQ  13) 

(EQ  14) 
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Treecode 


Figure  4.  Comparison  of  treecode  (top)  and  PIC  (bottom)  models  for  a 
one-dimensional  virtual  cathode  oscillation. 

Omelchenko  [8]  develops  a  3d-hybrid  PIC  code  using  the  assumptions  of  quasineutrality  and 
neglecting  displacement  current.  The  electrons  are  modeled  as  a  massless  fluid  collisionally  cou¬ 
pled  with  ions  and  neutrals.  The  ions  are  advanced  as  PIC  particles  with  additional  drag  terms  for 
collision  frequencies. 


V  x  B  =  pj 


(EQ  15) 
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Ve  +  V/-*- 


Je*B 

er]ec 


VPe 
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(EQ  16) 


FTa  “-y— ' v-M-pS-'.  +  Q.  <E<5|7> 

r(e,  are  the  electron  and  ion  resistivities,  respectively.  Qe  is  the  electron  heat  source,  and  y  is  the 
adiabatic  constant. 

Jones  et  al.  [9]  present  an  electromagnetic  PIC  program  utilizing  massless  electrons  and  including 
short-range  Coulomb  collisions.  For  interspecies  collisions,  the  fluid  properties  ( n ,  v,  T)  are  inter¬ 
polated  for  each  species  to  the  mesh  and  evaluated  as  a  collisional  force  field  when  integrating 
particle  motion  in  a  manner  that  satisfies  energy  and  momentum  conservation.  For  intraspecies 
collisions,  the  Langevin  equation  is  used  to  relax  particles  toward  a  Maxwellian  distribution.  The 
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collision  field  increases  computational  work  by  50%  and  allows  the  use  of  large  time  steps  rela¬ 
tive  to  the  collision  frequency.  The  electrons  are  modeled  as  massless  particles  with  the  assump¬ 
tion  of  quasineutrality: 


E  = 


-jexB-v(Pere) 

Pe 


m„ 


EZv«(v«“v<> 


(EQ  18) 


This  equation  is  solved  for  E  assuming  no  displacement  current.  The  equations  are  subcycled  to 
satisfy  the  vacuum  Courant  condition. 

In  addition  to  the  approximation  of  a  particle  species  either  partially  or  entirely  as  a  fluid,  another 
method  of  speeding  up  electromagnetic  simulation  is  using  an  implicit  solution  for  the  fields.  In 
an  explicit  electromagnetic  solution,  the  time  step  is  restricted  by  the  Courant  condition: 


A  t<—  (EQ  19) 

c 

This  is  often  considerably  more  prohibitive  than  the  stability  limitation  on  time  step  from  the 
explicit  particle  push: 


7 

At  <  —  (EQ  20) 

% 

While  not  a  fluid  method,  implicit  methods  must  often  be  coupled  with  fluid  models  to  obtain  sig¬ 
nificant  speed-up  and  can  also  provide  increased  computational  efficiency  to  particle  codes  with¬ 
out  fluid  species. 

Brackbill  and  Forslund  [10]  developed  an  implicit  discretization  of  Maxwell’s  equations  for  non- 
relativistic  plasma  simulation,  including  an  expression  for  current  using  the  plasma  pressure  ten¬ 
sor.  The  stability  is  limited  by  the  propagation  of  sound  waves  instead  of  electromagnetic  waves. 
High  frequency  waves  (those  not  resolved  by  the  time  step)  are  damped.  The  dispersion  in  the 
implicit  moment  equations  reduces  the  finite  grid  instability.  Accuracy  and  stability  are  empiri¬ 
cally  observed  when  At  satisfies 

(O(10_1)  <  vthAt/Ax  =  (\D/Ax)(®pAt))  <  0(1)  (EQ21) 

The  evolution  of  a  Weibel  instability  is  simulated  using  a  time  step  a>pAt  =  5.  No  energy  equation 
is  used  in  this  model. 

Bowers  [11]  developed  and  tested  an  implicit  electromagnetic  solver  appropriate  to  the  simulation 
of  high  frequency,  small  wavelength  plasma  systems.  An  implicit  Crank-Nicholson  finite  differ¬ 
encing  of  Maxwell’s  equations  is  used.  High  frequency  radiation  is  damped  and  performance  of 
the  field  solve  is  comparable  to  the  standard  explicit  solve.  Only  cases  with  modest  Courant  num¬ 
bers  retain  favorable  accuracy,  though  stability  is  retained  for  high  Courant  numbers.  The  method 
was  used  to  compute  small  wavenumber  dispersion  relations  of  bounded  plasmas. 

Relativistic  Monte  Carlo  Collision  Model 

Monte  Carlo  collisions  use  the  cross  section  to  calculate  the  probability  of  collision  and  the  scat¬ 
tering  angle.  The  method  of  calculating  collision  probability  is  the  same  as  for  non-relativistic 
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collisions,  so  the  null  collision  method  can  be  applied.  In  this  section,  we  consider  the  scattering 
angle  for  relativistic  collisions.  The  doubly  differential  cross  section,  Eq.  (22),  is  used  to  calculate 
the  scattering  angle  after  a  collision.  We  use  the  singly  differential  cross  section  to  get  the  second¬ 
ary  energy  after  an  ionization  collision. 

The  cross  section  for  the  relativistic  case  is  difficult  to  find  in  the  literature.  The  cross  sections  for 
the  non-relativistic  case  often  fail  for  the  intermediate  and  the  relativistic  regions.  Generally  the 
cross  sections  for  the  relativistic  case  do  not  work  well  for  the  other  cases.  We  define  the  criteria 
for  relativity  as  (3  =  vie  >  1/3  to  distinguish  the  collision  type  somewhat  arbitrarily;  one  could 
choose  a  different  boundary  when  more  precision  is  needed. 

For  relativistic  elastic  collisions,  we  will  use  the  Wentzel  Model  [12],  which  is  good  for  relativis¬ 
tic  collisions  and  also  is  reasonable  for  the  intermediate  regime.  The  angular  distribution  function 
of  the  Wentzel  model  is  as  follows: 


/(X) 


1  A(1  +  A) 

n(2A  +  1  -  cosx)2 


A.  =  i02(O.885Z  W3a0)  2(1.13  +  3.76(aZ/p)), 

P  (EQ  22) 

where  h  is  Plank’s  constant,/?  is  momentum,  aQ  is  Bohr  radius,  Z is  atomic  number,  and  a  =1/137 
is  the  fine-structure  constant. 

The  cumulative  distribution  for  this  distribution  function  gives  the  scattering  angle,  x 

cosx  -  2A  +  1  +  (EQ  23) 

where  0<  R<  1  is  a  uniformly  distributed  random  number. 

Using  the  basic  concept  for  the  relativistic  two  body  kinematics  and  the  scattering  angle  from  Eq. 
(22)  and  Eq.  (23),  the  energy  after  the  collision  is  obtained. 

For  the  elastic  collision  case,  we  need  the  general  equation  for  the  differential  cross  section.  After 
an  ionization  collision,  a  secondary  electron  is  created,  and  its  energy  must  be  computed.  The  sec¬ 
ondary  electron  energy  can  be  obtained  from  the  singly  differential  cross  section. 

The  calculation  method  is  the  same  as  for  the  non-relativistic  collision.  Only  the  equation  for  the 
singly  differential  cross  section  is  different.  The  equation  for  the  ionization  collision  cross  section 
is  [13]: 


*  (T,W,I) 


°(T)  f  1  |  4/  ,  4/  +  1 

f(T,I)\(W+P)2  3 (W  +  I)3  3(T-  W)3  ( T-W )2 


mc2(2T  +  me2) 


+ 


(W+I)(T-W)  2  2  2  2. 

\  /  ( T+mc  )  ( T+mc  ) 


) 


(EQ  24) 
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where  T  is  kinetic  energy  of  a  incident  electron,  W  is  secondary  electron  energy,  I  is  ionization 
threshold  and  is  normalizing  factor. 

We  use  the  cumulative  distribution  function  to  get  the  secondary  electron  energy  in  the  non-rela- 
tivistic  case. 


R  = 


(EQ  25) 


Once  we  get  the  secondary  electron  energy,  we  can  get  the  energy  of  two  electrons  and  one  posi¬ 
tive  ion.  Using  that  energy,  we  made  the  energy  conservation  equation. 


^ scat  ^ej  ^» 


Sinc  +  EN ' 


(EQ  26) 


where  escat,  s  and  zinc  are  energies  of  the  scattered,  ejected,  and  incident  electrons,  respectively. 
et  and  sA,  are  the  energies  of  the  created  ion  and  the  target  neutral  atom,  and  zion  is  the  ionization 
threshold  energy. 


Because  of  the  large  ion-to-electron  mass  ratio,  we  can  assume  that  the  momentum  of  the  incident 
electron  is  much  less  than  the  momentum  of  the  neutral  atoms.  In  other  words,  the  incident  elec¬ 
tron  strips  an  electron  off  the  neutral,  and  the  neutral  becomes  an  ion,  continuing  on  its  trajectory 
virtually  undisturbed.  This  assumption  allows  us  to  rewrite  the  above  energy  balance  equation  as: 


® scat  ^ej 


C  —  C  . 

inc  ion 


(EQ  27) 


Because  we  know  the  ionization  threshold  energy,  incident  energy,  and  ejected  energy,  we  can  get 
the  scattered  electron  energy  easily.  Once  we  get  the  scattered  electron  energy,  we  can  obtain  the 
scattering  angle  of  the  secondary  electron  in  a  same  way.  Obtaining  the  scattering  angle  requires 
the  doubly  differential  cross  section  as  in  other  collisions.  The  doubly  differential  cross  section 
that  is  used  for  the  relativistic  ionization  collision  is  [13]: 


<y(T,  W,%,I)cc  (, G3(T ;  W,  I)2  +  (cosx  -  G2(T,  W,  I))2)~\ 


where  G2(T,  W,  I) 


\(W+I)(T+2mc2) 
A|  T(W+I+2mc2)  ’ 


G3(T,W,I) 


a 


I(T -  W - 1) 
TW  ’ 


and 


a 


0.6 


(  2  \ 
\T  +  me  ) 


2 


(EQ  28) 


Next,  using  cumulative  distribution  function,  we  can  get  the  scattering  angle  after  the  ionization 
collision  event: 
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cos%  =  Gj  +  03  tan  (^(1  —  R)  atan  ^ —  R  atan 


l+G- 


(EQ  29) 


In  Phase  I,  the  relativistic  collision  model  was  implemented  in  the  ID  OOPD1  code. 


Neutral  Desorption  Model 

Partially  under  the  support  of  Phase  I,  a  model  for  the  surface  desorption  of  neutrals  was  devel¬ 
oped.  The  main  mechanism  for  desorption  due  to  energetic  ion  impact  appears  to  be  due  to  elec¬ 
tronic  sputtering,  in  which  energetic  ions  deposit  energy  in  a  medium,  creating  a  cylinder  of 
energetic  atoms  and  molecules  [14,15,16],  The  non-radiative  energy  relaxation  processes  trans¬ 
port  the  energy  back  to  the  surface  of  the  material,  where  it  leads  to  the  ejection  of  atoms  and  mol¬ 
ecules  when  the  binding  energy  is  exceeded.  Recent  molecular  dynamics  (MD)  results  have  been 
more  successful  than  so-called  “thermal  spike”  diffusive  models  for  these  sub-picosecond  phe¬ 
nomena. 


Consider  ion  impact  on  a  medium  at  an  angle  0  as  measured  from  the  surface  normal.  The  yield  at 
normal  incidence,  T0,  must  be  determined  empirically.  The  angular  dependence  is  obtained  from  a 
curve  fit  based  on  backscatter  data  from  [17]: 


^  =  1  +  1. 82  exp  (0.09x6).  (EQ30) 

*  0 

The  curve  fit  and  data  are  shown  in  Fig.  5.  The  backscatter  data  appears  to  correlate  well  with  the 
electronic  sputtering  rate,  rather  than  the  expected  1/cosO  dependence  from  energy  deposition 
considerations  [15]. 


The  energy  dependence  of  the  ejec¬ 
tion  distribution  is  given  by  [15]: 

9  TJF 

m  =  Cj - ^—3  (EQ  31) 

E,„(E+  Uf 

where  U  is  the  binding  energy, 

Eexc  is  the  excitation  energy,  and 
C[  is  a  normalizing  constant.  The 
energy  distribution  of  Eq.  (31)  is 
shown  in  Fig.  6.  This  is  imple¬ 
mented  in  the  model  by  inverting 
the  cumulative  distribution  func¬ 
tion  to  obtain  the  relation 


E 

U 


R  +  R 
1  +R 


(EQ  32) 


where  R  is  a  uniformly  distributed 
random  number,  0  <  R  <  1 . 


6 (deg) 


Figure  5.  Backscatter  coefficient  for  ion  impact 
based  on  data  from  [17]. 
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The  angular  distribution  of  the 
ejected  particles  is  chosen  based  on 
the  distribution  [15]: 

/(0)  =  C2cos(0)2,  (EQ  33) 

where  C2  is  a  normalizing  constant. 

The  angular  distribution,  which  is 
independent  of  the  incident  angle  as 
a  consequence  of  the  energy  relax¬ 
ation  mechanisms,  is  shown  in  Fig. 
7.  The  angle  is  chosen  from  this  dis¬ 
tribution  using  a  Monte  Carlo  rejec¬ 
tion  scheme. 

Task  2:  Electron  Scattering  from 
Metals/Insulators 


Figure  6.  Energy  distribution  of  atoms/molecules  for 
In  this  task,  the  focus  is  on  develop-  a  range  of  excitation  energies, 

ment  of  models  for  electron  scatter¬ 
ing  from  metals  and  insulators.  The  XOOPIC  code  currently  incorporates  the  Vaughan  model  for 
secondary  electron  emission  [18,19],  which  includes  both  single  event  backscattering  and  multi¬ 
ple  scattering  models.  Improvements  to  the  Vaughan  model  are  described  as  well. 


The  secondary  emission  coefficient  due  to  electron  impact  can  be  written: 


8(6.0)  =  Sm»0(i+*s6e2/27t  )»■(») 


(EQ34) 


Here,  the  incident  energy  is  given  by  8  ,  the 
angle  with  respect  to  the  surface  normal  is  0, 
ks§  is  a  surface  smoothness  parameter 
described  below,  k  is  also  described  below, 
and  8max  0  is  the  peak  coefficient  at  normal 

incidence  at  the  energy  8max  0 .  The  energy 

dependence  appears  implicitly  in  the  right 
hand  side  of  Eq.  (34)  through 


W(w)  =  j 


(wexp(l  -w))k,  w  <  3.6 
1.125/w0'35  ,  w  >  3.6 


,(EQ  35) 


where  the  normalized  energy,  w,  is  given  by:  Figure  7.  Angular  distribution  for  ejected 

particles. 

8—  8 

W  =  - °— - ,  (EQ  36) 

4axo(1+^n>e  /2<> 
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where  SQ  is  the  secondary  emission  threshold,  and  ksw  is  a  surface-smoothness  parameter  similar 
to  ks 5.  Both  ks §  and  ksw  vary  between  0  for  very  rough  surfaces  and  2  for  polished  surfaces.  Typi¬ 
cal  values  are  close  to  1.  The  exponent  k  in  Eq.  (35)  is  obtained  from  a  curve  fit  [19]: 


0.56,  w  <  1 
0.25,  1  <  w  <  3.6 


(EQ  37) 


The  energy  dependence  is 
shown  in  Fig.  8,  and  the  angular 

dependence  is  shown  in  Fig.  9  S  '  -v. 

for  electron  generated  secondar-  / 

ies.  Similar  yield  relations  can  0  / 

be  obtained  for  ion-induced  sec-  g  .  / 

ondaries,  sputtering,  and  des-  ^0.5  -  / 

orption.  •  / 

A  schematic  of  the  energy  dis- 
tribution  of  emitted  particles  is 

shown  in  Fig.  10.  Reflected  par-  [[,,,, . .  ,  ,  ,  . 

tides  typically  comprise  about  012345 

3%  of  the  total  particles  and  are  E/Emax0 

incident  electrons  reflected  at  Figure  8.  Energy  dependence  of  the  secondary  emission 
the  surface,  so  they  have  the  full  coefficient, 
incident  energy.  Backscattered 

primaries,  comprising  about  7%  of  the  ejected  electrons,  are  incident  electrons  that  have  pene¬ 
trated  the  surface  and  scattered  several  times  before  scattering  out  of  the  surface,  and  conse¬ 
quently  they  have  energy  between  zero  and  the  incident  energy.  True  secondaries,  comprising 
about  90%  of  the  ejected  population,  gain  sufficient  energy  to  escape  the  surface  by  interactions 
between  the  incident  particle,  the  free  and  bound  electron  populations,  and  the  solid  state  lattice. 
The  distribution  of  the  true  secondaries  is  taken  as  a  Maxwell-Boltzmann  distribution  at  a  temper¬ 
ature  Te  on  the  order  of  tens  of  eV: 
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/(■!)  =  - ^exp(-j-Q,<EQ38) 

where  is  the  Boltzmann  constant. 
The  angular  distribution  is  taken  to  be 
isotropic: 

g(0)  =  ~  cos(0) ,  (EQ  39) 

with  0  the  angle  from  the  surface  nor¬ 
mal.  The  azimuthal  angle  is  chosen 
randomly. 

In  the  Vaughan  model  for  secondary 
emission  used  in  XOOPIC,  a  step 
wise  value  of  k  is  used  (0.62  for  w  <  1 
or  0.25  for  w  >  1).  However,  another 
option  provides  an  energy-dependent 
value  for  A:  [19]: 


-90  -60  -30  0  30  60  90 


6  (degrees) 

Figure  9.  Angular  dependence  of  the  secondary 
emission  coefficient. 


k  = 


k^  k^ 
2 


k ,  -  k~, 


n 


atan(7rlnw) . 


(EQ  40) 


As  shown  in  Fig.  11,  when  an  energy-depen¬ 
dent  k  is  used,  based  on  Eq.  (40),  the  curve 
fits  the  experimental  data  better,  especially 
in  the  low  energy  regime  E  <  250  eV.  In  the 
high  energy  regime,  both  models  overesti-  f 
mate  the  secondary  yield. 


Figure  10.  Schematic  of  the  energy  distribution 
of  electrons  ejected  by  electron-surface  impact. 
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For  the  Vaughan  model,  the 
secondary  electron  emission 
coefficient  is  small  in  the  rela¬ 
tivistic  regime,  as  shown  in 
Fig.  12.  This  is  consistent  with 
experimental  data  in  the  litera-  | 
ture  for  the  MeV  range  [20],  | 

c 

In  Phase  I  of  this  work,  the  | 

XOOPIC  secondary  model  has  1 

been  improved  with  additional  | 

parameterization  so  that  the  | 

user  can  specify  the  full  set  of  § 
parameters  described  above 
from  the  input  file.  Further¬ 
more,  the  model  has  been  fur¬ 
ther  abstracted  so  that 
secondary  emission  is  now  a 
property  of  surfaces.  Under 

the  new  object  model,  each  Figure  1 1 .  Normalized  secondary  emission  coefficient  versus 

surface  may  have  multiple  impact  energy  (eV)  for  a  carbon  surface.  Experimental  data 

secondary  emission  properties,  from  Ritz  [28]. 
one  for  each  specie  of  incident 

particle,  with  independently  specified  energy  and  angular  dependent  yields  and  emission  distribu¬ 
tion  profiles. 


Impact  Energy 


Figure  12.  Secondary  emission  yield  in  the  relativistic 
regime  for  the  Vaughan  model. 
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In  a  previous  SBIR  program 
(DE-F G03  -00ER82964),  CCR 
funded  research  atN.C.  State 
University  to  develop  and  imple¬ 
ment  an  algorithm  to  accurately 
and  efficiently  predict  backscat- 
ter  yield  and  the  trajectories  and 
energies  of  electrons  backscat- 
tered  by  solids  [21].  The  algo¬ 
rithm  considers  the  energy  and 
direction  of  an  incident  electron, 
the  atomic  number,  atomic  mass, 
and  density  of  the  solid,  then 
determines  a  statistically  reason¬ 
able  path  for  the  electron 
through  the  solid  via  Monte 
Carlo  techniques.  This  so-called 
plural  scattering  model  traces  the 
path  of  an  electron  through  the 
solid  until  the  electron  loses  all 
of  its  energy  or  leaves  the  solid 
as  a  backscattered  electron.  Such  a  model  has  been  used  in  a  variety  of  applications,  but  in  this 
case  the  interest  is  predicting  the  behavior  of  backscattered  electrons.  When  applied  to  large  num¬ 
bers  of  electrons,  the  program  provides  statistically  accurate  results.  In  particular,  excellent  agree¬ 
ment  is  seen  between  the  backscatter  coefficients  measured  by  Hunger  and  Kuchler  and  those 
predicted  in  this  program  [22],  as  shown  in  Figure  13.  Furthermore,  the  angular  distribution  and 
energy  distributions  of  backscattered  electrons  predicted  are  consistent  with  those  measured  by 
Bishop  [23]. 


15  20  25  30 

Incident  Energy  (keV) 


35 


40 


45 


Figure  13.  Variation  of  the  backscatter  coefficient  with 
incident  energy  for  a  copper  target.  Each  predicted  value 
of  the  backscatter  coefficient  was  determined  by 
modeling  1000  electrons. 


In  Phase  I,  the  above 
model  was  implemented 
in  CCR's  Beam  Optics 
Analysis  (BOA)  finite 
element  code  to  track 
backscattered  electrons  in 
a  collector.  Figure  14 
shows  a  3D  conical  col¬ 
lector  with  three  genera¬ 
tions  of  secondary 
electron  trajectories.  Each 
generation  is  plotted  in  a 
different  color.  The  spent 
beam  is  injected  with 
appropriate  initial  condi¬ 
tions  in  an  external  mag¬ 
netic  field. 


Figure  14.  Simulation  of  a  collector  with  backscattered  electron 
emission.  Each  color  represents  one  generation  of  backscattered 
electrons 
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With  mesh  refinement,  the  mesh  density  is  increased  in  regions  containing  electrons  and  coars¬ 
ened  elsewhere.  Fig.  15  shows  the  final  mesh. 


Figure  15.  Unstructured,  adapted  mesh  in  a  collector  with  backscattered 
electron  emission. 


Task  3:  High  Voltage  Breakdown  on  Insulator  Surfaces 

Breakdown  due  to  high  voltage  at  insulator  surfaces  is  an  important  limiting  factor  in  achieving 
the  high  fields  needed  for  high  power  in  microwave  sources  [24].  The  effect  can  be  attributed  to  a 
number  of  interesting  physical  effects.  Here  we  consider  the  single  surface  multipactor,  possibly 
initiated  by  field  emission  (described  in  Task  1)  and  particle-impact  induced  outgassing  (desorp¬ 
tion  processes  described  in  Task  1). 

The  XOOPIC  code  includes  the  Vaughan  secondary  model  described  in  Task  2.  This  model  plays 
a  key  role  in  insulator  breakdown  by  single-surface  multipactor  discharge.  Single  surface  multi¬ 
pactor  has  been  described  in  the  literature  [25]  and  simulated  using  XOOPIC  [26].  In  the  single 
surface  multipactor  (see  Fig.  16),  electrons  in  the  vicinity  of  an  insulating  surface  are  accelerated 
in  RF  fields  in  the  direction  tangential  to  the  surface  normal,  gaining  energy.  A  net  positive  charge 
on  the  surface  supplies  a  DC  electric  field  component,  which  attracts  these  electrons  toward  the 
surface.  The  electrons  impact  the  surface  in  a  number  of  RF  cycles  determined  by  the  strength  of 
the  DC  electric  field  (and  hence  the  surface  charge  state),  the  initial  energy  of  ejection  for  the 
electron,  and  the  RF  frequency.  Upon  impact  with  the  insulator  surface,  the  electrons  may  gener¬ 
ate  secondary  electrons,  depending  on  the  energy  and  angle  of  impact.  Ejection  of  secondary  elec¬ 
trons  leaves  behind  an  immobile  positive  charge  at  the  insulator  surface,  which  generates  the  DC 
electric  field  that  causes  the  ejected  electron  to  make  an  excursion  in  the  direction  of  the  surface 
normal  and  eventually  return  to  the  surface.  When  the  gain  in  energy  during  an  excursion  due  to 
the  RF  fields  is  sufficient  to  result  in  impact  energies  and  angles  conducive  to  a  secondary  emis¬ 
sion  coefficient  greater  than  unity,  the  multipactor  can  grow  in  time.  Saturation  ultimately  occurs 
due  to  effects  at  the  edges  of  the  surface,  or  due  to  space  charge  effects  if  a  plasma  is  formed  at  the 
surface.  High  density  multipactor  can  result  in  absorption  and  reflection  of  RF  power,  which  can 
significantly  degrade  device  performance  or  even  damage  devices.  This  effect  can  occur  on  either 
side  of  the  surface. 
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Figure  16.  Schematic  of  insulator  breakdown  via  single  surface  multipactor. 


In  Phase  I  of  this  work,  the  field  emission  model  was  implemented  in  OOPD1  as  discussed  above, 
the  secondary  model  was  implemented,  and  the  Monte  Carlo  collision  model  was  extended  to  the 
relativistic  regime  as  described  in  Taskl  and  implemented  and  tested  in  OOPD1 . 

Also,  in  implementing  the  collision  model  in  OOPD1  and  verifying  with  XOOPIC,  a  discrepancy 
was  found  in  the  ion  drift  velocity  for  common  gases  with  the  same  cross  sections.  To  find  the 
drift  velocity  for  He  ion  in  helium  gas,  PIC/MCC  simulations  on  He  ions  were  performed  by 
using  XPDP1  and  XOOPIC  codes.  Two  parallel  electrodes  having  the  area  of  1  m2  are  separated 
by  the  gap  of  1  m.  The  voltage  of  200  V  is  applied  between  two  electrodes.  Initially,  the  ions  with 
uniform  density  and  zero  temperature  are  distributed  randomly  on  the  simulation  domain.  The  ion 
density  is  set  to  be  very  low  so  that  they  themselves  do  not  change  the  electric  field  of  200  V/m. 
The  pressure  of  background  helium  gas  is  1  Torr.  The  ions  gain  momentum  from  the  electric  field 
while  they  lose  momentum  through  collisions  with  background  gas.  After  several  collisions,  their 
macroscopic  mean  velocity  reaches  steady  state.  When  the  thermal  diffusion  of  ions  is  ignored, 
the  ion  drift  velocity  is  expressed  as 


(EQ  41) 


where  Mis  the  ion  mass  and  vm  is  the  ion-neutral  collision  frequency  for  momentum  transfer.  In 
XPDP1  and  XOOPIC  codes,  the  elastic  and  charge  exchange  collisions  for  ions  are  taken  into 
account. 


Fig.  17  shows  the  ion  drift  velocity  obtained  from  recent  XPDP1  and  XOOPIC  simulations.  Since 
the  simulations  were  performed  for  identical  physical  systems,  the  drift  velocities  should  be  the 
same,  but  differ  by  a  factor  of  8.  In  the  XOOPIC  simulations,  the  steady-state  ion  drift  velocity 
varied  for  different  ion  densities  or  number  of  superparticles. 
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From  these  indications,  we  cor¬ 
rected  errors  in  the  Monte  Carlo 
source  code  of  XPDP1  and 
XOOPIC.  To  validate  the  modi¬ 
fied  codes,  the  theoretical 
velocity  of  Eq.  (41)  was  com¬ 
pared  with  that  obtained  from 
XPDP1  and  XOOPIC  for  an 
analytically  tractable  model 
elastic  collision  cross  section 
that  is  inversely  proportional  to 
the  ion  velocity  (so  that  the  col¬ 
lision  frequency  is  independent 
of  the  ion  velocity).  For  sim¬ 
plicity,  the  charge  exchange 
collision  was  neglected.  The 
simulation  results  of  XPDP1 
and  XOOPIC  are  now  in  good 
agreement  with  the  theoretical 
result,  as  shown  in  Fig.  18. 


Figure  17.  Erroneous  result:  the  drift  velocity  as  a  function 
of  time  does  not  agree. 


Simulations  using  experimental 
cross  sections  also  show  good 
agreement  between  XPDP1  and 
XOOPIC  simulation  results  for 
identical  cross  section,  without 
dependence  on  density  or  numer¬ 
ical  particle  weight,  as  shown  in 
Fig.  19. 

Task  4:  Intense  Heat  Fluxes 

In  high  voltage  microwave 
sources,  intense  heat  fluxes  occur 
when  energetic  electrons  impact 
surfaces  such  as  collectors, 
anodes,  or  beam  intercept  at 
unintended  locations.  Other 
mechanisms  of  heat  deposition 
can  include  plasma  bombard¬ 
ment.  Although  electron  flux  ide¬ 
ally  should  be  confined  to  the  collector  surface,  electrons  often  also  impact  other  surfaces.  Beam 
instabilities  can  deposit  significant  fluxes  on  unintended  surfaces,  including  the  anode  and  the 
slow  wave  circuit,  resulting  in  heating  that  leads  to  outgassing  or  even  permanent  damage.  The 
outgassing  can  result  in  plasma  formation,  which  further  disrupts  the  beam  and  absorbs  or  reflects 
microwave  power,  decreasing  device  efficiency  or  disrupting  high  power  operation. 


Time  (us) 

Figure  18.  Good  agreement  for  modified  XPDP1  and 
XOOPIC  with  theoretical  result  for  analytic  model  cross 
section  with  energy-independent  collision  frequency. 
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Since  the  heat  flux  corresponds 
to  energy  deposition  by  parti¬ 
cles,  a  kinetic  method  of  track¬ 
ing  particles  and  measuring  the 
energy  deposited  at  surfaces  is 
required.  Since  sputtering  and 
surface  damage  are  also  of 
interest,  it  is  useful  to  model 
the  distribution  of  particle 
impacts  in  space,  time,  energy 
of  impact,  and  angle  of  impact. 

XOOPIC  currently  includes  the 
capability  to  measure  the 
energy  and  angle  of  impact  for 
a  given  surface  separately,  in 
addition  to  measuring  the  cur¬ 
rent  collected  on  surfaces  for 
each  species  in  the  simulation. 

The  existing  model  is  currently 

limited  to  simple  single  element  surfaces  and  provides  neither  time  resolution  nor  spatial  resolu¬ 
tion.  In  particular,  the  code  currently  collects  the  following  de-correlated  quantities  for  impacts: 


Figure  19.  Modified  XPDP1  and  XOOPIC  show  good 
agreement  of  ion  drift  velocity  for  experimental  cross 
sections. 


In  Phase  I,  a  model  for  obtaining  the  fully  correlated  surface  flux,  /(x,0,£),  was  investigated.  The 
model  will  allow  the  user  to  select  some  subset  of  the  fully  correlated  flux  distribution  for  activa¬ 
tion  on  individual  surfaces,  in  order  to  minimize  computation  and  memory  requirements. 


The  impacting  distribution  function  of  particle  species  is  known  explicitly  in  the  PIC-MCC  meth¬ 
odology.  The  impacting  energy,  E,  is  known  from  the  particle  velocity,  and  the  angle  from  a  sur¬ 
face  normal,  0  is  known: 


0  =  atan(— )  (EQ42) 

Vvll 

The  particle  energy  and  angle  may  be  weighted  to  a  mesh  when  impacting  boundaries  of  interest 
to  obtain  a  gridded  representation  of  the  impacting  distribution  function,  J(E,Q).  Fluxes  to  the  sur¬ 
face  (energy  flux,  momentum  flux,  etc.)  may  be  calculated  from  the  distribution  function.  The 
average  power  absorbed  by  the  surface  over  time  period  T  is 

\  j5  [0/2  f;°xdEdMSf(E,  0)  (EQ  43) 

The  integral  \  dS  is  over  surface  area.  In  ID  codes,  the  integration  reduces  to  a  multiplication.  For 
2D  and  3D  codes, /is  additionally  a  function  of  position  on  the  surface  and  is  4-  and  5-dimen- 
sional,  respectively.  An  example  of  a  distribution  function,/  is  shown  in  Fig.  20. 

XOOPIC  collects/©)  and  flE)  independently.  This  saves  memory,  as  two  dimensions  of  data  are 
stored  per  boundary  (f  is  a  function  of  position  along  boundary  and  E  or  0)  compared  to  three 
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Figure  20.  Example  distribution  function, 

dimensions  for J(E,Q)  and  position.  However,  since/)))  and/(£)  are  uncorrelated,  the  fluxes  (heat 
flux,  momentum  flux,  etc.)  cannot  be  calculated.  Thus  it  is  desirable  to  obtain/is,0). 

The  memory  needed  to  store / depends  on  what  information  is  stored.  If  the  full,  time-dependent 
distribution  function  is  stored  in  bins  (f(x,t,E, 0)),  the  memory  needed  to  store  the  data  is 


d- 1 

<eq44> 

i 

In  this  equation,  Mis  the  memory  used,  Nx,i  is  the  number  of  grid  points  on  the  surface  in  the  i  - 
direction,  and  d  is  the  number  of  dimensions  in  configuration  space.  Nq  and  NE  are  the  number  of 
bins  for  angle  and  energy,  respectively,  and  NT  is  the  number  of  time-steps  stored.  Min  Equation 
(44)  must  be  multiplied  by  the  number  of  species,  though  this  is  usually  small.  Averaging/ over 
time  reduces  the  system  dimension,  reduces  the  needed  memory  by  a  factor  of  Np  and  reduces 

noise  as  an  affective  ensemble  average.  Another  way  to  reduce  the  needed  memory  is  to  use  a 
functional  (spline)  representation  of /instead  of  bins,  which  may  be  considered  zero-order 
splines.  Higher  order  splines  may  be  used  to  accumulate  a  smooth  distribution  function  that  has 
one  or  more  continuous  derivatives  and  conserves  one  or  more  moments  of  the  distribution  func¬ 
tion  (number,  momentum,  energy,  etc.).  Note  that  binning/ only  conserves  the  lowest  moment 
(number).  Using  a  spline  representation  of /  conserves  memory,  reduces  noise,  and  provides  a 
convenient  way  to  integrate  the  distribution  function  to  calculate  fluxes.  However,  the  coefficients 
of  the  spline  must  be  recalculated  as  particles  hit,  which  reduces  computational  efficiency.  The 
trade-offs  between  computational  efficiency  and  adequate  resolution  of  the  distribution  function 
will  be  studied.  The  binned  data  may  also  be  post-processed  to  reduce  noise  and  obtain  a  spline 
representation  of '/  but  higher  moments  of  the  distribution  function  are  not  conserved. 
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The  use  of  a  single  angle,  0,  assumes  azimuthal  symmetry  about  this  angle,  which  is  usually  ade¬ 
quate  to  represent  particle  interaction  with  an  isotropic  surface.  An  additional  angle,  VF,  may  be 
introduced  if  the  full  distribution  function  needs  to  be  calculated.  However,  this  increases  compu¬ 
tational  time  and  storage  requirements,  so  it  should  be  implemented  with  an  option  to  enable  the 
full  calculation  only  when  required. 

Task  5:  X-ray  Generation  and  Effects 

X-rays  cause  a  number  of  concerns  for  high  voltage  electron  devices.  First,  they  pose  a  safety  haz¬ 
ard  to  personnel  and  sensitive  equipment,  requiring  shielding  and  other  precautions.  X-rays  can 
cause  photoemission  from  various  surfaces  inside  microwave  devices,  leading  to  electrons  in 
undesired  areas,  which  can  interact  with  microwaves,  background  gas,  and  the  electron  beam, 
adversely  impacting  device  performance  and  possibly  even  disrupting  device  operation.  The  elec¬ 
trons  generated  by  photoemission  can  also  lead  to  dual  surface  and  single  surface  multipactors, 
which  can  cause  significant  reflection  and  absorption  of  microwave  power  as  well  as  space  charge 
effects.  X-rays  can  also  photoionize  background  gas  and  surface  impurities,  leading  to  plasma  for¬ 
mation,  even  far  from  the  electron  beam.  Indeed,  cross  sections  for  ionization  can  be  much  larger 
for  electrons  formed  in  modest  voltage  regions  than  for  high  voltage  beam  electrons,  so  that  these 
photoelectrons  can  become  a  more  significant  source  of  plasma  than  the  more  dense  beam. 

Energetic  charged  particles  undergoing  rapid  acceleration,  for  example  in  collisions  with  neutral 
particles,  interaction  with  strong  magnetic  fields,  or  collisions  with  solid  matter,  emit  electromag¬ 
netic  radiation  in  a  broad  spectrum,  including  x-radiation.  In  such  problems  it  is  useful  to  develop 
the  formalism  in  a  way  that  relates  the  radiation  intensity  and  polarization  directly  to  properties  of 
the  charged  particle  trajectory  and  motion.  The  total  radiation  emitted,  the  angular  distribution  of 
radiation,  and  its  frequency  spectrum  are  also  of  interest.  For  non-relativistic  motion  the  radiation 
is  described  by  the  Larmor  result,  and  for  relativistic  motion  we  have  more  sophisticated  theory 
and  calculations  [27]. 

The  Larmor  power  emission  for  a  non-relativistic,  accelerated  charge  is  as  follows: 


P  = 


(EQ  45) 


Larmor’s  formula  Eq.  (45)  can  be  converted  into  a  suggestive  form  to  find  the  appropriate  gener¬ 
alization. 


„  _  2_e _ ( dg 

~  3w2c3 \dt  '  dt. 


(EQ  46) 


where  m  is  the  mass  of  the  charged  particle,  p  its  momentum.  The  Lorentz  invariant  generaliza¬ 
tion  is 


P  = 


2  g2  (dPydp^ 

3  rnc'dx  dx  ) 


(EQ  47) 


where  dx  =  dth{  is  the  proper  time  element,  and  is  the  momentum-energy  four-vector. 
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The  following  is  the  four-vector  scalar  product: 


(EQ  48) 


dx  dx 


In  linear  acceleration,  the  motion  is  one-dimensional.  From  Eq.  (48)  it  is  evident  that  in  that  case 
the  radiated  power  is 


P 


2  e2  f dp\2 

3  2  AdxJ 
Jm  c  Ul 


(EQ  49) 


The  rate  of  change  of  momentum  is  equal  to  the  change  in  energy  of  the  particle  per  unit  distance. 
Consequently 


2  e2  fdE\2 

^  32  3V/&; 

Jm  c  ux 


(EQ  50) 


showing  that  for  linear  motion  the  power  radiated  depends  only  on  the  external  forces  that  deter¬ 
mine  the  rate  of  change  of  particle  energy  with  distance,  not  on  the  actual  energy  or  momentum  of 
the  particle. 

In  circular  accelerators  like  the  synchrotron  or  betatron,  circumstances  change  drastically.  In  such 
a  situation,  the  momentum  p  changes  rapidly  in  direction  as  the  particle  rotates,  but  the  change  in 
energy  per  revolution  is  small.  This  means  that: 


dj) 

dx 


,,  1  dE 

HpI  »  -- 


Then  the  radiated  power  Eq.  (47)  can  be  written  approximately 

P 


2  e2  22,  ,2  2 e2c  4n4 

y  co  |p!  =  3— y  (3 

Jm  c  p 


where  we  have  used  ©  =  c(3/p,  p  being  the  orbit  radius. 
The  radial  component  of  Poynting’s  vector  can  be  written: 


S  n  = 


2  f  1 
e  1 

4  7tc  p2 


n  x  [(n  -  (3)  x  g] 


(1-P-ny 


(EQ  51) 


(EQ  52) 


(EQ  53) 
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In  Eq.  (53)  S  •  n  is  the  energy  per  unit  area  per  unit  time  detected  at  an  observation  point  at  time  t 
'  of  radiation  by  the  charge  at  time  f  =  t-R{t')/c .  If  we  wanted  to  calculate  the  energy  radiated 
during  a  finite  period  of  acceleration,  from  f  =  Tx  to  f  =  T2,  we  would  write 


T2  +  [R(T2)/c] 
r,  +  [R{T,)/c] 


[S  ■  n ]dt  =  l 


=  t2 
-Ti 


(S  n  )£* 


We  therefore  define  the  power  radiated  per  unit  solid  angle  to  be 


(EQ  54) 


-  i?2(Sn)J;  =  *2S-n(l-P-n)  (EQ 55) 

If  we  assume  the  charge  to  be  accelerated  only  for  a  short  time  during  which  (3  and  (3  are  essen¬ 
tially  constant  in  direction  and  magnitude,  and  we  observe  the  radiation  far  enough  away  from  the 
charge  that  n  and  R  change  negligibly  during  the  acceleration  interval,  then  Eq.  (55)  is  propor¬ 
tional  to  the  angular  distribution  of  the  energy  radiated.  With  Eq.  (53)  for  the  Poynting  vector,  the 
angular  distribution  is 


dP(t')  =  e2  |n  x  [(n  -  ft)  x  (3"j| 2 
dQ.  4  nc  (i_p.„)5 


(EQ  56) 


If  0  is  the  angle  of  observation  measured  from  the  common  direction  of  (3  and  f3 ,  then  Eq.  (56) 
reduces  to 


dP(t’)  _  e2v2  sin20 
^  27rc3(l  -  Pcos0)5 

For  small  angles,  the  angular  distribution  Eq.  (57)  can  be  written  approximately 


(EQ  57) 


dP(t')  .  8  e2v2  8 
dQ  %  c3  7 


(Y0)2 

(i+y2e2)5 


(EQ  58) 


Another  example  of  angular  distribution  of  radiation  is  that  for  a  charge  in  instantaneously  circu¬ 
lar  motion  with  its  acceleration  (3  perpendicular  to  its  velocity  p.  We  choose  a  coordinate  system 

such  that  instantaneously  p  is  in  the  z  direction  and  p  is  in  the  x  direction.  With  the  customary 
polar  angle  0,  <j)  defining  the  direction  of  observation,  as  shown  in  Fig.  21,  the  general  formula  Eq. 
(56) reduces  to 


dP(f) 

dQ 


4tic3(1  -  Pcos0)3 


.2  2 
sin  0  cos  t|) 

y2(l  -  Pcos<j>)2 


(EQ  59) 


Using  the  above  equations,  we  can  treat  each  radiation  as  a  photon  with  a  discrete  energy  [12]. 
The  radiation  is  negligible  beyond  the  critical  frequency, 


Virtual  Laboratory  Environment  for  High  Voltage  Radiation  Source  ExperimentsMay  2,  2005 


25 


(EQ  60) 


®erf/  =  3T3f- 
rL 

where  rL  is  the  Larmor  radius. 

Consider  a  1  GeV  electron  orbiting  a  10  T  static  magnetic 
field.  The  relativistic  momentum  factor  is  y  =1952,  and  the 
velocity  is  effectively  the  speed  of  light,  c,  to  six  places.  The 

gyroradius  is  0.334  m.  The  power  radiated  is  P  =  6x  10'6  W, 
so  the  electron  would  lose  1%  of  its  energy  in  about  270  ns. 

1  o 

The  critical  frequency  is  cocrit/27r  =  3.2x  10  Hz,  with  a  cor¬ 
responding  critical  wavelength  of  X,crit  =  9.4x  10"11  m.  In 
order  to  resolve  the  wavelength  on  the  mesh,  we  would  need 

Ax<  A,/4  =  2.4x  10’ 11  m.  For  a  nominal  system  length  of  L  = 

3rL  in  each  dimension,  the  required  number  of  cells  in  each 

dimension  is  4.2  x  1010,  or  1.7  x  1021  cells  in  a  two  dimen¬ 
sional  simulation.  This  exceeds  the  memory  capacity  of  exist¬ 
ing  computers.  Furthermore,  the  time  step  required  to  satisfy  Figure  21 .  Coordinates  for 

,  relativistic  collision  dynamics, 

the  Courant  condition  is  A  t  <  A  xN2c  =  5.7  x  10'zu  s.  The 

cyclotron  period  is  7  ns,  so  one  orbit  would  require  1.2x  1011 

time  steps,  while  running  to  1%  energy  loss  would  require  4.7x  10  time  steps.  Resolving  the 
critical  frequency  requires  a  time  step  A  t  <  3.1x  10'19s,  which  is  even  more  untenable. 

The  conclusion  of  this  analysis  is  that  the  radiation  damping  cannot  be  self-consistently  modeled 
by  the  standard  PIC  scheme  even  for  a  single  particle.  Due  to  the  grid  resolution  and  time  step 
constraints,  there  is  no  possibility  of  directly  simulating  the  fields  of  the  radiation  emitted  over  a 
spectrum  sufficient  to  include  a  significant  fraction  of  the  radiated  power. 

However,  from  a  previous  AFOSR-funded  STTR  collaboration  between  U.C.  Berkeley  and 
TechXCorp.,  a  model  was  implemented  in  XOOPIC  to  account  for  the  energy  lost  by  particles  due 
to  emission  of  radiation,  using  an  equivalent  radiation  damping  force: 


where  p  is  the  particle  momentum,  and 


(EQ  61) 


(EQ  62) 


This  radiation  damping  term  allows  for  the  power  loss  of  accelerating  particles,  even  though  it  is 
not  possible  to  track  the  radiation  emitted  on  these  time  and  space  scales. 
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Task  6:  Detailed  Structural  Visualization 

The  VLAB  code  will  be  based  on  an  extended  version  of  UCB’s  XOOPIC  code.  This  code  is  a 
very  flexible  and  powerful  analysis  tool  for  particle  simulations,  but  it  is  difficult  for  new  and 
casual  users,  because  it  lacks  a  graphical  user  interface  (GUI).The  problem  geometry  and  run  time 
information  is  specified  using  a  text  file  of  the  required  format.  There  are  several  hundred  input 
parameters  documented  in  forty  pages  of  text.  This  presents  a  daunting  learning  curve  for  new  or 
casual  users.  In  the  Phase  II  program,  we  propose  to  develop  a  GUI  that  will  simplify  the  program 
for  the  new  and  casual  user  while  not  hindering  the  advanced  user. 

The  proposed  GUI  will  significantly  ease  the  generation  of  the  data  required  for  XOOPIC.  The 
primary  improvement  will  come  from  graphical  display  and  input  of  the  problem  geometry.  The 
problem  parameters  will  be  grouped  logically  into  the  categories  of  Setup,  Boundaries,  Ports, 
Emitters  and  Diagnostics.  Advanced  features  of  OOPIC,  such  as  Variables,  will  be  supported  and 
enhanced  with  user-defined  groups.  Finally  a  powerful  calculator  based  diagnostics  capability 
will  be  added  to  significantly  increase  the  data  analysis  capabilities.  These  features  and  others  will 
be  described  in  the  following  sections. 

The  GUI  will  be  developed  using  the  Qt  application  development  framework.  The  Qt  application 
program  interface  (API)  is  consistent  across  Windows,  Linux,  Unix,  and  Macintosh  operating  sys¬ 
tems,  thus  enabling  a  platform  independent  application  deployment.  Initial  releases  of  the  GUI 
will  target  the  Windows  and  Linux  environments. 

Main  Screen 

The  start  up  screen  will  be  similar  to  the  one  shown  in  Figure  22.  If  an  input  file  is  selected,  a 
graphical  output  of  the  geometry  will  be  displayed  in  the  problem  geometry  window.  During  a 
computation,  output  of  particle  positions  can  be  also  flagged  for  display  on  this  screen.  The  five 
main  menus  provide  a  logical  grouping  of  the  many  potential  parameters  controlling  the  simula¬ 
tion.  They  are  grouped  as  follows: 

•  Setup 

•  Variables  -  user  defined  symbolic  variables 

•  Grid  -  dimensions  of  simulation  and  grid  details 

•  Control  -  time  step  size,  field  solver  type 

•  Species  -  particle  definitions 

•  Description  -  text  description  of  simulation 

•  Boundaries 

•  Conductor  -  perfect  conductor  boundaries 

•  Equipotential  -  perfect  conductor  boundaries  that  are  at  a  constant  potential 

•  Dielectric  -  dielectric  boundaries 

•  Polarizer  -  a  boundary  that  allows  fractional  particle  transmission 

•  Current  -  region  of  current  flow 

•  Ports 

•  Exit  -  a  region  with  minimal  reflection  of  electromagnetic  fields 

•  Gap  -  a  field  excitation  boundary 

•  Gaussian  -  a  Gaussian  profile  electromagnetic  pulse  excitation  port 
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•  Emitters 

•  Beam  -  region  to 
generate  particles 
defined  by  species 
types 

•  Variable  weight  - 
spatial  variable 
weight  emitter 
beam  emitter 

•  Fowler-Nordheim  - 
field  emission  sur¬ 
face 

•  Particle  Creation 

•  Collisions  -  Monte 
Carlo  collision 
parameters 

•  Load  -  spatial  and 
Maxwellian  veloc¬ 
ity  distributions  of 
named  particle  spe¬ 
cies 

•  Plasma  source  -  rectangular  region  for  plasma  generation 

•  Diagnostics 

•  Default  -  list  of  default  diagnostics  for  display 

•  Calculator  -  user  defined  operations  on  analysis  variables  for  display 

Geometry  Entry 

In  the  current  version  of  XOOPIC,  the  geometry  of  the  simulation  is  entered  as  a  series  of  mesh 
coordinates  (or  physical  coordinates  that  are  mapped  to  the  nearest  mesh  coordinates).  The  text 
entry  of  these  coordinates  is  very  laborious,  prone  to  error,  and  provides  no  immediate  visual 
feedback  for  error  checking.  In  the  proposed  GUI,  the  input  can  be  entered  using  the  current 
method  of  mesh  or  physical  coordinates  along  with  a  new  option  of  entry  using  cursor  position  on 
the  mesh  display.  An  example  of  coordinate  entry  is  shown  in  Figure  23.  Selecting  a  point  on  the 
screen  will  enter  the  grid  and  physical  coordinates  for  the  grid  point  nearest  the  selected  point.  If 
the  user  enters  the  grid  point  by  typing  in  the  grid  point  box,  the  physical  coordinate  will  be 
entered  automatically  in  the  corresponding  box. 

In  the  entry  example  depicted  in  Figure  23,  the  input  is  for  entry  of  a  single,  linear,  line  segment. 
For  generation  of  more  complex  geometries,  drawing  element  types  such  as  polyline,  arcs,  parab¬ 
olas  and  circles  will  be  supported.  These  drawing  primitives  will  allow  a  user  to  quickly  define 
the  problem  boundaries  for  the  majority  of  simulation  types;  however,  we  will  also  investigate 
construction  of  a  translator  to  utilize  geometries  generated  in  2D  CAD  programs.  This  translator 
would  allow  complex  geometries  to  be  drawn  using  the  advanced  drawing  facilities  available  in  a 
CAD  package  familiar  to  the  user.  This  drawing  can  then  be  exported  in  industry  standard  formats 
such  as  ACIS  for  import  to  the  proposed  XOOPIC  GUI.  CCR  implemented  this  capability  into  its 
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Figure  22.  Proposed  main  screen  for  OOPIC  GUI. 
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Figure  23.  Entry  of  coordinate  information  for  a  new  boundary. 


current  Beam  Optics  Analysis  product.  After  the  CAD  drawing  is  imported  by  the  GUI,  the 
boundary  types  will  be  identified  for  generation  of  the  simulation  parameters. 

Most  boundary  types  require  entry  of  a 
only  a  few  parameters  to  specify  the 
boundary  conditions.  In  addition  to 
boundary  specific  parameters,  however, 
there  are  a  large  number  of  generic 
boundary  conditions  that  can  be  speci¬ 
fied  but  are  typically  left  at  their  default 
values  or  do  not  apply  to  the  current 
boundary  type.  The  GUI  will  separate 
the  required  entries  from  the  generic 
entries  using  tables  selected  by  tab  but¬ 
tons.  In  addition,  only  those  generic 
boundary  conditions  that  are  applica¬ 
ble  to  the  current  boundary  will  be  dis¬ 
played.  An  example  boundary  dialog 
screen  is  shown  in  Figure  24.  This 
screen  also  shows  how  the  particle  spe¬ 
cies  will  be  selected  from  the  list  of 
defined  types.  The  species  types  will  be 
a  database  of  predefined  types,  such  as  electrons,  and  the  GUI  will  allow  entry  of  other  user 
defined  types. 


Figure  24.  Example  dialog  box  for  entry  of 
secondary  boundary  condition. 
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Variables 


XOOPIC  allows  user  defined  symbolic 
variables.  This  facility  is  used  to  create 
definitions  that  can  be  referenced  as  input 
for  other  simulation  parameters.  The  GUI 
will  support  this  feature  and  also  allow  the 
variables  to  be  collected  in  named  groups. 
An  example  dialog  screen  showing  the 
Variable  implementation  is  shown  in 
Figure  25.  In  this  example,  three  groups  of 
variables  are  defined  (Constants,  Grid 
Variables  and  Plasma  Variables)  that  are 
accessible  by  selection  of  the  appropriate 
tab  button.  The  GUI  will  allow  creation, 
modification,  and  deletion  of  each  variable 
group  and  variables  in  the  group. 


Constants  Grid  Variables  Plasma  Variables 


plasmaDensity  [2.1  e+20 

simVolume  jpi  *  maxRadiusMKS  “  maxRadiusIvIKS 
numRcIsPerCell  Js 

thermalSpeed  |  spe  edOfLi  g  M  *  sqrt(  p  lasiriaDen  sily  )| 


Edit  Variables 


Edit  Groups  OK 


Diagnostics 


Figure  25.  Variable  definition  dialog  box. 


The  XOOPIC  code  graphically  displays  a  large  number  of  default  diagnostics  that  are  updated  as 
the  simulation  progresses.  Typical  diagnostics  include  particle  position  or  velocity  vs.  coordinate, 
charge  density,  Poynting  calculation,  and  field  quantities.  In  addition,  user  defined  diagnostics  are 
available.  This  feature  allows  generation  of  time  histories  of  simulation  quantities  such  as  the 
field,  Poynting  vector  and  charge  density  along  a  line,  point  or  rectangular  spatial  region.  Addi¬ 
tional  operations  that  can  be  performed  on  analysis  data  include  integration,  polarization  calcula¬ 
tions  and  power  spectral  density  calculations. 


The  GUI  will  retain  all  the  existing  diagnostic  plotting  capabilities  of  XOOPIC  and  significantly 
enhance  the  user  defined  diagnostic  capability.  An  arc  element  will  be  added  to  the  geometries 
(line,  point,  rectangle)  that  can  be  used  for  analysis  quantities.  The  operations  that  can  be  per¬ 
formed  on  these  quantities  will  be  greatly  increased.  A  listing  of  the  operations  that  will  be 
included  are: 

•  General  operations:  absolute  value  and  smoothing, 

•  Complex  operations:  Re,Im  and  conjugate, 

•  Scaler  operations:  reciprocal,  power,  square  root,  trig  operations,  derivatives  and  log, 

•  Vector  operations:  dot  and  cross  products,  divergence,  curl,  tangent  and  normal. 


In  addition  to  generating  plots,  an  option  to  export  to  file  will  be  added. 


Summary 

In  summary,  the  Phase  I  program  achieved  all  the  initial  objectives  and  demonstrated  feasibility 
for  all  tasks. 
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